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The detection of gravitational waves from extreme-mass-ratio inspirals (EMRIs), comprising a 
stellar-mass compact object orbiting around a massive black hole, is one of the main targets for low- 
frequency gravitational-wave detectors in space, like the Laser Interferometer Space Antenna (LISA) 
or evolved LISA/New Gravitational Observatory (eLISA/NGO). The long-duration gravitational- 
waveforms emitted by such systems encode the structure of the strong field region of the massive 
black hole, in which the inspiral occurs. The detection and analysis of EMRIs will therefore allow 
us to study the geometry of massive black holes and determine whether their nature is as predicted 
by General Relativity and even to test whether General Relativity is the correct theory to describe 
the dynamics of these systems. To achieve this, EMRI modeling in alternative theories of gravity is 
required to describe the generation of gravitational waves. However, up to now, only a restricted class 
of theories has been investigated. In this paper, we explore to what extent EMRI observations with 
a space-based gravitational- wave observatory like LISA or eLISA/NGO might be able to distinguish 
between General Relativity and a particular modification of it, known as Dynamical Chern-Simons 
Modified Gravity. Our analysis is based on a parameter estimation study that uses approximate 
gravitational waveforms obtained via a radiative-adiabatic method. In this framework, the trajectory 
of the stellar object is modeled as a sequence of geodesies in the spacetime of the modified-gravity 
massive black hole. The evolution between geodesies is determined by flux formulae based on 
general relativistic post-Newtonian and black hole perturbation theory computations. Once the 
trajectory of the stellar compact object has been obtained, the waveforms are computed using the 
standard multipole formulae for gravitational radiation applied to this trajectory. Our analysis is 
restricted to a five-dimensional subspace of the EMRI configuration space, including a Chern-Simons 
parameter that controls the strength of gravitational deviations from General Relativity. We find 
that, if Dynamical Chern-Simons Modified Gravity is the correct theory, an observatory like LISA 
or even eLISA/NGO should be able to measure the Chern-Simons parameter with fractional errors 
below 5%. If General Relativity is the true theory, these observatories should put bounds on this 
parameter at the level £ 4 / 4 < 10 4 km, which is four orders of magnitude better than current Solar 
System bounds. 

PACS numbers: 04.30.Db, 04.30.-w, 04.50.Kd, 95.30.Sf, 97.10.Sj 



I. INTRODUCTION 

There is strong observational evidence for the exis- 
tence of black holes in galactic x-ray binary systems, 
seen as ultraluminous x-ray sources, and in the centers 
of galaxies, seen as active galactic nuclei (see, e.g. [1]). 
Indeed, observations carried out by space- and ground- 
based telescopes suggest the presence of a dark compact 
object, likely a massive black hole (MBH), at the center 
of most observed galaxies (see [2] and references therein). 
In a typical galaxy, the MBH is surrounded by around 
10 7 — 10 s stars forming a cusp or core (see, e.g. [3]). As 
a consequence of relaxation, mass segregation and large 
scattering encounters between the stars, stellar compact 
objects (SCOs) may be perturbed onto orbits that pass 
sufficiently close to the MBH and become gravitationally 
bound forming a binary system. Therefore, the capture 
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of a SCO by a MBH is likely to be a frequent phenomenon 
in the Universe. 

Once the SCO has become bound to the MBH, it starts 
a slow inspiral driven by the emission of gravitational 
waves (GWs). During this process, the system loses en- 
ergy and angular momentum and the orbit of the SCO 
circularizes and shrinks adiabatically, i.e. on a timescale 
much longer than the orbital period. The loss of en- 
ergy and angular momentum occurs initially in bursts, 
when the object passes through the orbital pericenter, 
but eventually the gravitational radiation is being emit- 
ted continuously until the object reaches the innermost 
stable orbit and plunges into the MBH. For EMRIs whose 
GW frequencies lie in the sensitivity band of space-based 
GW detectors, like LISA (Laser Interferometer Space An- 
tenna HE]) or eLISA/NGO (evolved LISA/New Gravi- 
tational Observatory [fOE]), the central MBH must have 
mass in the range, M # ~ 10 4 — 10 7 M Q . The systems 
of interest must also have a SCO compact enough to 
avoid tidal-disruption and so the SCO must be a stel- 
lar mass black hole (m + sw 1 — 50Mq), a neutron star 
(m 4 « 1.4M©) or a white dwarf (m 4 « O.6M ). The 
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typical mass ratios, fi = m^/M,, of EMRI systems are 
therefore in the range ~ 10~ 6 — 10~ 4 . 

The strongest detectable EMRI signals are unlikely to 
be any closer than a luminosity distance D ~ 1 Gpc [S] , 
at which distance the instantaneous amplitude of the 
measured EMRI signal is an order of magnitude below 
the level of instrumental noise and the GW foreground 
from galactic white-dwarf binaries. EMRI detection will 
therefore rely on matched filtering of the detected data 
stream with a bank of templates of the possible signals 
that might be present in the data. During the last year 
before plunge, an EMRI will generate ~ 1/fj, gravita- 
tional waveform cycles in the LISA band [9]. During 
this time, the orbit of the SCO tracks the strong field 
geometry in the vicinity of the MBH and maps out the 
(multipolar) structure of the MBH spacetime [10] in the 
emitted GWs. 

GWs from EMRIs are generated in the strong field re- 
gion close to the MBH and therefore probe General Rel- 
ativity (GR) in a regime which, up to now, has not been 
reached observationally. If GR is the true theory of grav- 
ity describing EMRI dynamics, their waveforms will de- 
termine the parameters of the system with very high pre- 
cision. However, if the central MBH is not described by 
the Kerr metric or GR does not properly describe the bi- 
nary dynamics in the strong field regime, and we assume 
GR when constructing our detection templates, we will 
obtain incorrect results from GW observations. There- 
fore, there is a strong motivation for studying what kind 
of modifications to the dynamics of EMRIs one could ex- 
pect from considering well-motivated theories of gravity 
other than GR. To that end, we must understand how 
the signals are modified in these alternative theories so 
that we are able to detect and quantify deviations from 
GR. 

The use of EMRI observations for such tests of fun- 
damental physics has been explored by several authors 
(see [11] and references therein), but the majority of that 
work has focussed on using the observations to constrain 
the properties of "bumpy" black holes. These are so- 
lutions to the field equations of general relativity that 
represent spacetimes that differ from the Kerr solution 
by an amount controlled by a tunable deviation param- 
eter. EMRI observations will be able to place bounds 
on the size of deviations of the forms considered [I2HI6] . 
However, this is not necessarily a test of general relativ- 
ity, since the bumpy black holes are constructed within 
that theory. It is rather a test of the "no-hair" prop- 
erty of black holes (stationary astrophysical black holes 
are described by the 2-parameter (mass and spin) fam- 
ily of spacetime geometries of Kerr [17] ) and hence the 
auxiliary assumptions that go into the no-hair conjec- 
ture. Hence, "bumpy" black holes are actually a test of 
the Kerr geometry assuming GR is the correct theory of 
gravity. 

Due to the myriad of alternative theories of gravity 
available, the questions that arise are: Which kind of 
theory do we choose to compare against? What new fea- 



tures might we expect to observe in the GW signals that 
might allow us to distinguish this theory from GR? In 
this paper we address these different questions and ex- 
plore the capability of a space-based detector like LISA 
to discriminate between GR and an alternative theory 
of gravity. In particular, we focus in a modification 
of GR constructed by the addition of a Chern-Simons 
(CS) gravitational term (also known as the Pontryagin 
invariant) to the action. Interest in this theory was ini- 
tiated with the work of Jackiw and Pi |18| where grav- 
itational parity violation was investigated. Such a term 
appears in four-dimensional compactifications of pertur- 
bative string theory due to the Green-Schwarz anomaly- 
canceling mechanism [19j and also in loop quantum grav- 
ity when the Barbero-Immirzi parameter is promoted to 
a scalar field coupled to the Nieh-Yan invariant [2"0Tl22) . 
Moreover, the Pontryagin term is unavoidable in an effec- 
tive field theory (see [23] in the context of cosmological 
inflation). In the approach of Jackiw and Pi, the Pon- 
tryagin term is introduced in the action multiplied by a 
scalar function and, in this way, it contributes to the field 
equations (in a four-dimensional spacetime the Pontrya- 
gin term is a topological invariant and hence does not 
contribute to the field equations), but this field is not 
dynamical. That is, it is a given function of the space- 
times coordinates. This version of CS modified gravity 
has been extensively studied and it has been shown to 
be dynamically too restrictive and, for instance, generic 
oscillations of non-rotating Schwarzschild black hole are 
not allowed [Mj- In addition there are problems with 
the uniqueness of solutions of the theory [25] . For these 
reasons we focus on the version of the theory in which 
the CS scalar field is dynamical, i.e. Dynamical Chern- 
Simons Modified Gravity (DCSMG), see [5B] for a review 
of CS modified gravity. 

The first study of EMRIs in DCSMG was done in [27], 
where the main ingredients of the problem were discussed 
and a simple waveform model was put forward. This 
model used the so-called semi-relativistic approximation, 
in which the trajectories are geodesies and the waveforms 
are built by using a standard multipolar expansion of 
the gravitational radiation. Then, differences between 
the GR and DCSMG waveforms were studied and also 
some predictions for the relative dephasing of the waves 
were made. However, this work relied on the assump- 
tion that radiation reaction (RR) effects, i.e. the effects 
that arise from the interaction of the SCO with its own 
gravitational field would allow one to distinguish between 
GR and DCSMG. Without RR the harmonic structure of 
the waveforms is going to be very similar and hence it is 
likely that it would be always possible to match a signal 
with both GR and DCSMG template waveform models. 
On the other hand, in a recent study [2S], corrections to 
the gravitational- and scalar-wave fluxes for circular or- 
bits around a non-rotating MBH in CS gravity have been 
computed using perturbation theory. This type of com- 
putations are very promising and can complement the 
work we present in this paper. 
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In this paper we go beyond the model of [37J by in- 
cluding two important additional ingredients: (i) RR ef- 
fects based on a hybrid scheme [53] that combines (post- 
Newtonian) PN approximations and fits to Teukolsky re- 
sults |30| ; (ii) Fisher parameter estimation techniques to 
make predictions on the capability of a space-based de- 
tector to measure the EMRI parameters, in particular 
a CS parameter that controls the deviations from GR. 
We have built kludge waveforms in the spirit of [3T] and 
have used them to estimate expected measurement er- 
rors for the main parameters describing an EMRI system 
in DCSMG. We find that for LISA these error estima- 
tions have the following order of magnitude: central black 
hole mass, A log M. ~ 5 • ; central black hole spin, 
Aa ~ 5 ■ 10~ 6 M. ; orbital eccentricity, Ae ~ 3 ■ 1CT 7 ; 
luminosity distance of source, Alog(Z?i//i) ~ 2 • 1CV 2 ; 
and for the CS parameter, £, in the combination ( = a£, 
we find A log £ ~ 4 • 10~ 2 . Moreover, we also use this 
framework to put bounds on the CS parameter, £, di- 
rectly. Assuming that GR is the correct theory to de- 
scribe ERMIs, we find that LISA measurements could 
put bounds of the order t; 1 / 4 < 10 4 km, which are bet- 
ter by four orders of magnitude than those derived from 
frame dragging observations around the Earth |32) . 

This paper is organized as follows. In Section [TT] we 
describe all the components used for the construction of 
EMRI gravitational waveforms in DCSMG and the re- 
sponse of space-based GW detectors. This includes the 
basic aspects of the theory, the deviations in the MBH 
geometry and its impact in the orbital dynamics and the 
inclusion of RR effects. In Section ITTT1 we summarize the 
basics elements of signal analysis theory and parameter 
estimation based on Fisher matrix techniques. In Sec- 



tion IV we apply these techniques to the waveforms and 
response models built in Sec. [TTJ providing parameter er- 
ror estimates for both LISA and eLISA/NGO and also 
bounds to the CS parameter. We finish in Section VI 
with conclusions and a discussion. Appendix [A] con- 
tains the form of the power spectral density of LISA and 
eLISA/NGO, while Appendix [B] contains the formulae 
needed for the construction of the RR effects. 

Throughout this paper we use Einstein summation 
convention for repeated indices and geometrized units 
in which G = c = 1. Spacetime indices are denoted 
by Greek letters; spatial indices are denoted with Latin 
letters . . .; denotes the canonical metric covari- 
ant derivative operator and □ = g^V^V,, denotes the 
d'Alambertian wave operator. 



II. EMRIS IN DCSMG 

In order to carry out parameter estimation studies to 
assess the ability of a given GW detector to detect and 
extract the physical information of an EMRI system, we 
first need a theoretical model of the generated waveforms. 
EMRIs are complex systems and we do not have yet a 
description accurate enough to produce waveforms in GR 



that can be used for data analysis purposes. However, 
for parameter estimation studies it is enough to have a 
waveform model that contains all the features of the real 
waveforms and that approximates the waveform phase to 
within a few cycles over the whole inspiral. 

Due to the large difference between the masses of the 
two components in an EMRI, the GW signal can be mod- 
eled accurately using perturbation theory (see e.g. |33|1. 
where the SCO is represented as a structureless particle 
orbiting in the MBH spacetime background. Although on 
short timescales the orbit of the SCO is approximately 
a geodesic of the MBH spacetime, its parameters slowly 
change with time due to RR effects. The best method 
we have to estimate these RR effects is the so-called 
self-force approach. At present, the gravitational self- 
force has been computed for the case of a non-rotating 
MBH [34j [35] and progress is being made towards cal- 
culations for the more astrophysically relevant case of a 
spinning MBH [3S] (see [3TH5U] for reviews). 

In parallel to the self-force program, some efforts to 
build certain approximation schemes to model EMRIs 
have been made. For the purposes of this work we focus 
on the so-called Numerical Kludge waveform model |31) . 
In that framework, the orbital motion is given by a se- 
quence of geodesies around a Kerr MBH, with the evo- 
lution of the geodesic parameters dictated by a dissipa- 
tive RR prescription. This prescription is based on PN 
evolution equations for the orbital elements (from 2PN 
expressions for the fluxes of energy and angular momen- 
tum) calibrated to more accurate Teukolsky fluxes with 
45 fitting parameters [2S] . The waveforms are then mod- 
eled using a multipolar expansion [4"U] . 

To accurately compute the GW emission from EMRIs 
in an alternative theory of gravity, we need to understand 
both how the orbital dynamics of the binary are altered 
and how gravitational wave generation and propagation 
differs in the alternative theory. In DCSMG, the GW 
emission formulae are not modified at leading order [27] , 
and so in this paper we will consider modifications to the 
underlying orbital dynamics only. In what follows we de- 
scribe the main components of our waveform model, sum- 
marizing the procedure introduced in |27j and including 
the RR effects just described. 



A. Formulation of DCSMG 

In DCSMG the action functional depends on the space- 
time metric g^,,, on the CS scalar field i?, and on the 
matter fields ip mat , and it can be cast in the following 
form 

+ 2 S *lE»u>0] + Smatfe/^! VVatl ) (!) 

where k n is the gravitational constant, l/(167r) in ge- 
ometrized units, and a and j3 are universal coupling con- 
stants that control the strength of the CS modifications. 
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The different contributions to the action are: 
Einstein-Hilbert action 



^eh — 



the GR 



(2) 



where g is the metric determinant and R is the Ricci 
curvature scalar; the CS gravitational correction 



stress-energy tensor has the same form (in terms of the 
gauge- invariant metric perturbation describing the GWs) 
as the one of Isaacson for GR. This is due to the fact that 
the averaging involved in the short-wave approximation 
cancels out all the CS corrections giving rise, at leading 
order, to essentially the same backreaction in DCSMG as 
in GR. 



^cs 



RR , 



(3) 



where * RR := *R C 



l5 RP„z = he^"R a a ,.,.R P ^ is the 



Pontryagin density, R^ va a is the Riemann tensor, 
is the Levi-Civita antisymmetric tensor and here the as- 
terisk denotes the dual operation; the CS scalar field ac- 
tion term 

' (V„0) (V„0) + 2V(i?)] , (4) 

where V is the scalar field potential, which is neglected 
in this work (i.e. V = 0); and finally, S mat [g^,tp mat } is 
the action of the different matter fields. 

Varying the action with respect to the metric and the 
CS scalar field, we obtain the field equations of DCSMG: 



G, 



K 



a = J_ 



rpxi 



a 



find = -^*RR, 
4 



(5) 
(6) 



where G^ v is the Einstein tensor and C^ v is the so-called 



C-tensor which has two parts, C^ v = Cf" - 
Cf = {VJ)e aSv ^ a V v R p \, 



C%" with 



°2 



(7) 



Finally, T^f* is the matter stress-energy tensor and T^ v 
is the stress-energy of the CS scalar field, given by 



t(#) 







(V^)(V^) 
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One can see that taking the divergence of the field equa- 
tions ([5]), using the Bianchi identities and the conserva- 
tion of the matter stress-energy tensor, one obtains the 
field equation (|6| for the CS scalar field. 

There are several consequences of DCSMG that are rel- 
evant for this work. The first one is that the number of 
independent waveform polarizations that a detector far 
away from a GW source will see are the same in DCSMG 
as in GR [57], i.e., the plus and cross tensor GW polariza- 
tions. In addition to the plus and cross polarizations, in 
DCSMG there is an additional breathing mode, however 
it decays faster, typically like r~ 2 and is therefore unlikely 
to be detected by an observer far away from the source. 
Another important property of GWs in DCSMG is the 
structure of the stress-energy (or mass) tensor that can 
be associated with the GWs in the short-wave approx- 
imation (see, e.g. [H]), commonly known as the Isaac- 
son tensor g3] . In [27] it was shown that the GW 



B. The MBH geometry in DCSMG 

The first ingredient we need to model the dynamics of 
an EMRI system is the geometry of the MBH. In GR we 
know that, provided the no-hair conjecture is true, all 
MBHs must be described by the Kerr metric. However, 
this is no longer true in DCSMG. We do not have an exact 
solution in DCSMG for spinning MBHs, but there is an 
approximate solution [25 , 44J that has been found using a 
small-coupling approximation (using £ cs = a 2 /(M # /3«: N ) 
as the expansion parameter, with M. being the MBH 
mass) and a slow-rotation approximation (defined by 
a/M, < 1, with a = \S m \/M„ < a/M, < 1, and S. 
is MBH spin). Using a system of coordinates that in the 
GR limit coincide with the well-known Boyer-Lindquist 
(BL) coordinates (t, r, 9, 4>) [5S], the non- vanishing metric 
components have the following form 



1 



2M.r 



p_ 
A ' 



sm 



it</> 



5 £ a M. 

8 Mf W.~r~^ 



12M. 27 Ml 



Tr- 



lOr 2 



2M.ar 



sin 2 6» , 



(9) 

(10) 
(11) 
(12) 



(13) 



where we have introduced the following definitions: p 2 — 
r 2 + a 2 cos 2 6, A = r 2 f + a 2 , f = 1 - 2M./r, and 
S = (r 2 + a 2 ) 2 -a 2 Asin 2 #. The effects of the CS gravita- 
tional modification are parametrized by a single universal 
constant, £, given by 



a 

0kI 



(14) 



Notice that the only metric component that gets mod- 
ified with respect to the general relativistic case is the 
component g t( j [Eq (13)]. The term in this component 
that is proportional to the CS parameter £ falls off with 
distance as r -4 , that is, it decays much faster than the 
rest of the metric components and hence it becomes neg- 
ligible at large distances. Only gravitational systems like 
EMRIs can probe this modification as they penetrate into 
the strong field region of the MBH. 



5 



At the level of approximation at which the DCSMG 



metric [Eq ( 13 ) ] was obtained, it is possible to show that 



it has most of the properties of the Kerr metric [27] , 
in particular the DCSMG metric is stationary and ax- 
isymmetric, and also has a Killing tensor, which is im- 
portant for an analysis of the orbital motion. Moreover, 
the DCSM metric has the same algebraic structure as 
the Kerr one. Regarding the multipolar structure of 
this DCSMG metric, let us remember that the multi- 
pole moments of the Kerr metric are fully determined by 
the MBH mass and spin (or equivalently, by the mass 
monopole and current dipole) according to the following 
simple relations: Me+iSg — M (ia) e , where {Mg} e=0 x 
and {Sg}g =Q x are the mass and current multipole mo- 
ments respectively. The multipole moments associated 
with the CS metric deviate from those of Kerr starting 
at the S 4 multipole, as one can see by employing the 
multipolar formalism of [20] (see also [55]). Despite this 
deviation, the structure of these multipole moments still 
preserves the philosophy of the no-hair conjecture since 
they only depend on the mass and spin of the MBH. 
There is also a dependence on the CS parameter £, but 
this is a universal dependence that would be the same for 
all MBHs and hence it cannot be considered to be hair 
of the MBH. 

The equations for the metric and CS scalar field are 
coupled, so they have to be solved simultaneously. The 
solution for the CS scalar, at the same level of approxi- 
mation as for the metric, is: 



■6 



5 a a cos ( 



8 M. 



1 



2M. ISMl 



5r 2 



(15) 



This scalar field falls off as r 
finite energy associated with it. 



and therefore it has a 



C. Orbital Kinematics 

It was argued in [27. that in DCSMG massive parti- 
cles should follow geodesies of the spacetime metric. At 
the lowest order of approximation, and for short periods 
of time, the trajectory of the SCO can then be approxi- 
mated by geodesies of the metric given in Eqs. 
Actually, we are going to approximate the orbital motion 
as a sequence of geodesies, as in the GR case within the 
Numerical Kludge (NK) waveform model. For this rea- 
son, it is important to analyze in detail the structure of 
the geodesic motion around this modified-Kerr metric. 

In the previous subsection we mentioned that the mod- 
ified MBH geometry has essentially the same physical and 
geometrical properties as the Kerr metric. In particular 
it has the same number of symmetries. Therefore, we can 
separate the geodesic equations as in the Kerr case, intro- 
ducing certain constants of the motion. More specifically, 
we have the energy per unit SCO mass, E, the angular 
momentum component along the spin axis per unit SCO 
mass, L z , and finally the Carter constant per unit SCO 
mass squared, C . 



The geodesic equations have the following struc- 
ture E7J: 



i = i K + L z 5g^ s (r), 
i> = i> K -E8gf{r), 



2EL z f5g$ s (r), 



9 2 = 2 K , 



(16) 
(17) 
(18) 
(19) 



where the dots denote differentiation with respect to 
proper time. The quantities (t K , f K , 9 K , <p K ) are the 
counterparts of the geodesic equations in the Kerr met- 
ric, which are given by (see, e.g. |46j ) 



p 2 i K = —a(aE sin 2 — L x ) 

m 2 i „2 

[(r 2 + a 2 ) E - aL z ] 



A 



(20) 



P^K = ~ [(r 2 + a 2 )E - aL z ] (aE ^ (21) 

p 4 r 2 K =[(r 2 +a 2 )E-aL z } 2 

~A[Q+{aE-L z ) 2 +r 2 ], (22) 

P 4 ^k = Q~ cot2 eL l ~ a 2 cos2 0(1 - E<1 ) » ( 23 ) 

where Q is an alternative definition of the Carter con- 
stant, related to C by 



= C-(L Z - aE)" 



(24) 



It is clear from Eqs. (H6H18) that the CS deviations 



are determined by a single function of the radial coordi- 
nate r, Sg9 s (r), which has the form 



7? S = ^ 

H 112r 6 / 



70 + 120- 



189- 



Mi 



(25) 



Only the equation for the polar coordinate is un- 
changed. Since we are dealing with bound orbits, both 
the radial and the polar motion include turning points 
(extrema of motion) at which the time derivatives, r or 
6, vanish. This can create numerical problems when 
integrating the set of Ordinary Differential Equations 
(ODEs) given by Eqs. ( 16 1-( 18 1. To avoid this, we follow 
the same strategy as in the case of Kerr geodesies and in- 
troduce two angle coordinates, ip and \: associated with 
the radial and polar motion respectively: 



P M, 



1 + e cos tp 



cos 2 9 = cos 2 6> min cos 2 x , (26) 



where p and e are the dimensionless semilatus rectum 
and the eccentricity of the orbit respectively, and # min is 
the minimum of 9 in the orbit (the turning point in the 
polar motion). The orbital parameters (e,p) are related 
with the radial turning points, the apocenter (r apo ) and 
pericenter (r ,), through the standard expressions: 



pM, 



pM, 



(27) 



6 



or equivalently 



D. Orbital Dynamics 



2r 



M.(r 



(28) 



The radial coordinate r oscillates in the interval 
(r pcri , r apo ). Similarly, given the turning point of the po- 
lar motion, 6* min e [0,7r/2], performs a liberation mo- 
tion in the interval (0 min ,7r — # mil J). We introduce the 
orbital inclination angle (with respect to the spin direc- 
tion), through the following relation 



0i„ c = si & n ( L z) 



(29) 



where sign(L z ) = 1 corresponds to a prograde orbit and 
sign(L 2 ) = — 1 corresponds to a retrograde orbit. A dif- 
ferent definition of the orbital inclination angle can be 
given in terms of the constants of motion (E, L z , C or Q) 



cos i 



Q 



(30) 



In general, both inclination angles, 6 lac and l, are quite 
similar |47) and coincide in the non-spinning limit, a = 0. 

We work with two geodesic parametcrizations, one 
based on the orbital parameters [e,p,9 inc or t) and one 
based on the constants of motion (E, L z , C or Q). Chang- 
ing from one parameterization to the other is a funda- 
mental step in our computations. In the GR case there is 
a well-known procedure (see, e.g. [33 SB] ) to do so. Here, 
we have used the implementation described in the appen- 
dices of However, these formulae are only valid in 
GR and, in our case, the CS modification of the radial 
equation of motion changes the location of the turning 
points. In practice this translates into a different rela- 
tion between the two sets of parameters (e,p, # inc or t) 
and (E, L z , C or Q). Given that we are not dealing with 
large deviations from the GR case, we have used a numer- 
ical procedure based on the Newton-Raphson method for 
finding roots (see, e.g. 50 ), where the values obtained 
from the GR method have been used as the starting point 
for the iteration algorithm. We have seen that in practice 
this works quite well and the iteration converges rapidly 
to the correct values. 

Finally, due to the separability of the geodesic equa- 
tions, which is closely related to the spacetime symme- 
tries, we can distinguish in the motion three fundamen- 
tal frequencies (here with respect to coordinate time t) 
associated with the radial motion, f r = l/T r (T r is the 
average time to go from pericenter to apocenter and back 
to pericenter), with the polar motion, f e — 1/Tg (T s is 
the average time for a full oscillation of the orbital plane, 
going from 9 — min to 9 = tt — 9 mia and back to 9 — 6* min ), 
and = 1/T, (T, is the average time for the SCO's az- 
imuthal angular coordinate <fi to cover 2tt radians). It is 
important to mention that these frequencies change when 
including the CS modifications [27] ■ 



So far, the orbital dynamics described have been for 
geodesic orbits. In order to compute the SCO trajectory, 
we gradually evolve the parameters of the instantaneous 
geodesic orbit under RR. The RR effects not only drive 
the SCO inspiral, but also break the degeneracy between 
orbits in GR and others in DCSMG since a given ini- 
tial orbital configuration will evolve differently in these 
theories. Therefore we need to implement RR effects in 
the EMRI dynamics in the framework of DCSMG. 

As we have mentioned above, in this paper we adapt 
the NK waveform model to the case of DCSMG [S]. In 
the NK waveform model, the RR driven evolution uses a 
"hybrid" scheme described in 29 , where formulae for the 
evolution of the constants of motion (E, L Z ,C or Q) are 
derived in terms of PN approximations (at 2PN order) 
combined with fits to results from the Teukolsky formal- 
ism (see [301 151)). In principle one should then derive 
the analogous formulae for the case of DCSMG, but this 
involves a number of major developments that are cur- 
rently out of reach. Instead, we will take into account 
one of the important results about GWs in DCSMG dis- 
cussed previously — the realization that the stress-energy 
momentum tensor for GWs, the Isaacson tensor, has the 
same form in terms of the GW metric perturbation in 
both theories, GR and DCSMG. This means that, to 
leading-order, the properties of the GW emission in GR 
and DCSMG are the same. Then, we approximate the 
fluxes of energy and angular momentum in the GWs, 
and also the evolution of the Carter constant under GW 
emission, by using the GR expressions. In what follows, 
we describe the formulae and procedures to update the 
geodesic orbits in our NK-EMRI model. 

The evolution equations for the constants of motion 
(E,L Z ,C or Q) have the following structure: 



dE 
~dt 

dt 

dQ 

dt 



= M/E( a >P. e >0tao) > 
= A 1 /L i ( a )-P> e ^i„c) , 

= tJ-fo( a ,P,e,9 iiic ) . 



(31) 
(32) 

(33) 



The evolution equation for C follows from these equa- 
tions and Eq. (24). The form of the right-hand sides 
f E , f L , and Jq of Eqs. ( [3lj)-(33) and full details of 
their derivation can be found in [29 . In Appendix [Tj] we 
summarize the main expressions needed to build these 



1 Without RR effects the phase evolution of an EMRI waveform, 
in both theories, will be a multiple Fourier series of the three 
fundamental frequencies, i.e. it will contain harmonics of the type 
exp{27ri/ m n p t} with f m>n>p = mf r +nf g +pf <fi . Then, we would 
be able to associate physical parameters with a given EMRI in 
both theories, and hence we would not be able to discriminate 
between them. 
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right-hand sides and thus evaluate the evolution of the 
constants of motion. 

In practice, there are two ways in which we can use 
the evolution equations (31)-(33). The first one consists 



of computing a phase-space trajectory for the orbital pa- 
rameters by integrating the set of ODEs for the evolu- 
tion of the energy, E, angular momentum component 
along the spin axis, L z , and Carter constant, Q. Once 
the phase space trajectory (E(t), L z (t), Q{t)) has been 
computed, these time-dependent constants are used on 
the right-hand side of the geodesic equations (20)- (23 1, 



to construct the inspiral trajectory of the SCO in the 
Boyer-Lindquist-like coordinates of the MBH spacetime. 
The second option, the one that we use in this paper, is 
to consider the extended system of ODEs consisting of 
Eqs. |20l)-((23l and Eqs. |3ll)-((33| and integrate them to- 



gether in time. Although this system of ODEs is coupled, 
there is a clear hierarchical structure, since the subsystem 
of Eqs. ( |31[ )-(33) can in principle be integrated indepen- 
dently of the subsystem of Eqs. ( 20 )-( 23 ) , which can be 
seen as a subsidiary system. 

As mentioned before, Eqs. (31)-(33) are in principle 



only valid in GR. If the true theory of gravity is DC- 
SMG, these evolution equations will contain corrections. 
At leading order, GW emission in DCSMG takes the 
same form as in GR [57], but corrections to the fluxes 
will still arise from the DCSMG modifications to the 
orbital motion. These corrections were computed for 
circular orbits in DCSMG in [52], but enter at a high 
post-Newtonian order. For this reason, we do not make 
any modifications to the GR expressions but directly em- 
ploy the fluxes described in Gair and Glampedakis [29] . 
Although we therefore use the same RR formulae to 
evolve the trajectory in DCSMG as in GR, this still leads 
to a different SCO evolution, since the dependence of 
the orbital elements (e,p, 6> inc or i) on the "constants" 
of motion (E,L Z ,C or Q) is different in the two the- 
ories, which leads to correspondingly different gravita- 
tional waveforms. That is, the mapping between the or- 
bital elements (e,p, 6 inc or t) and the constants of motion 
(E,L Z1 C or Q) is different in DCSMG and GR. There- 
fore, given some initial orbital parameters (e ,p , 9 inc or 
l ), after evolving the EMRI system for some time, the 
final orbital parameters in GR will be in general different 
from the orbital parameters in DCSMG. 

Taking into account the previous considerations, the 
inspiral is constructed in the following way: For a 
given set of initial orbital parameters (e o ,p o ,0 inc or 
l ), we find the associated initial constants of the mo- 



tion (E Q , L z ,C or Q ), which are different from the 
ones that we would obtain in GR for the same ini- 
tial eccentricity, semilatus rectum and inclination an- 
gle. Subsequently, we evolve the "constants" of motion, 
(E\ ,L Z \ ,Q\ ), under RR, using the method described 
above. Then, from the current values of the constant of 
motion, (E , L z , C or Q ), their rates of change due to 

RR, (E\ ,L Z \ ,Q\ ), and the value of the radial period, 
T r (the time to go from the apocenter to the pericenter 



and back again to apocenter) [2"T] 
constants of motion, (E 1 , L z 1 ,Q 1 ) 
equations: 



we obtain the new 
using the following 



E Q + E\ N r T r , 



Qi = Qo + Q\oKT r , 



L z \ N r T r , 



(34) 
(35) 
(36) 



and N r is a pre-specified parameter that represents the 
number of radial periods elapsed between each update of 
the constants of the motion. The expression for C 1 fol- 
lows from these formulae for (E 1 ,L z l ,Q 1 ) and Eq. (24|. 
Finally from (E 1 ,L Z 1 ,C 1 /Q 1 ), we obtain the new val- 
ues of the orbital parameters (e 1 ,p 1 ,6'. nc i/ix). This al- 
gorithm is iterated along the whole EMRI evolution to 
obtain the SCO orbit. In Figure[l]we illustrate a section 
of a generic orbit for a typical system that we use later 
in our parameter estimation analysis (see Table [n| . 




FIG. 1. Depiction of a section of the inspiral orbit for an 
EMRI (system A in Table [TTf with parameters M. = 5 • 
10 5 M o , a = 0.25M., e = 0.25 and ( = 5 ■ 10" 2 M. 5 . 



E. Waveform Modeling and Detector Responses 

In the previous subsections we have seen how the tra- 
jectory of the SCO is obtained and, in the following, we 
describe how we compute the gravitational waveforms 
and the response of the LISA and eLISA/NGO detec- 
tors. Following [3T] and [37] we employ the multipo- 
lar expansion of the metric perturbations describing the 
GWs emitted by an isolated system, which assumes that 
the GWs propagate in a flat background spacetime to 
reach the observer/detector |40j . In this work, we con- 
sider only the lowest-order term, the mass quadrupole. 
This term involves second time derivatives of the trajec- 
tory and these are readily obtained from the geodesic 
equations (20 1- (23 1. Then, the transverse-traceless (TT) 
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metric perturbation is computed from the following ex- 
pression 



(*) = Ih 



(37) 



where 1^ denotes the mass quadrupole and r the luminos- 
ity distance from the source to the observer. In terms of 
the source stress-energy tensor, T , the mass quadrupole 
moment is: 

n STF 



d A xx % x 3 T u (t,x') 



(38) 



where STF stands for symmetric and trace-free. Treat- 
ing the SCO in the point-mass approximation, the non- 
vanishing components of the stress energy tensor have 
the following form: T u (t,x l ) = pit.x 1 ) and T tj {t,x i ) = 
p{t,x l )vi {t), where p(t,x l ) is the energy density of the 
SCO which, in the point-mass limit, is given by: 

p{t,x l ) = mj {3) [**-*'(*)] , (39) 

where 5^ denotes the three-dimensional Dirac delta dis- 
tribution, z l (t) are the spatial Cartesian coordinates (as- 
sociated with the flat spacetime background) of the SCO 
trajectory and v l {t) — dz l (t)/dt are the components of 
the corresponding spatial velocity. To evaluate this in our 
model we make a "particle-on-a-string" approximation, 
that is we identify the Boyer-Lindquist-like coordinates 
(r, 9, <f)) of the SCO's orbit with flat-space spherical polar 
coordinates, and introduce Cartesian coordinates in the 
usual way 

x = r sin 9 cos (f> , y = r sin#sin0 , z = r cos f? .(40) 

Although this description leads to inconsistencies, like 
the non-conservation of the flat-space energy-momentum 
tensor of the particle motion, it has been found to work 
well when generating EMRI waveforms in GR [31] and 
so we do not expect this to introduce large errors in the 
waveforms, in particular in the phase. A possible alterna- 
tive could be to use coordinate systems more adapted to 
the multipolar expansion of the gravitational radiation, 
like harmonic or asymptotic-Cartesian mass-centered co- 
ordinates [3D] (see also [4"9"]). 

We now consider detection of these signals by a space- 
based detector in a heliocentric orbit (like LISA or 
eLISA/NGO). We describe the direction from the de- 
tector to the EMRI system by a unit 3-vector n, which 
also gives the propagation direction of the GWs from the 
EMRI to the detector. The orthogonal plane to n is the 
GW polarization plane and we can introduce there two 
unit and orthogonal vectors p and q by using the spin 
direction, S. = aM m z: 



P = 



n x z 
In x z| 



q=pxn. 



(41) 



The vectors (n, p, q) form a spatial orthonormal basis 
that can be used to construct the GW polarization ten- 
sors: 



PiPj - mj > 



(42) 



The corresponding plus, h + , and cross, h xi GW polar- 
izations are given by: 

1 



M*), (43) 



2 + xv ' 2 

and the complete GW metric perturbation is: 

M*) = 4M*) + e «M*)- ( 44 ) 



Using Eqs. (37 1 and (|39|) we obtain the following simpli- 



fied expressions for the GW polarizations in terms of the 
SCO position z t (t), velocity v l (t) — dz l /dt, and acceler- 
ation a l (t) = d 2 z l /dt 2 : 

2m+ 



:(*) 



+ .x 



[a i {t)ss i (t)+v i (t)v j (t)] ■ (45) 



Once we have the GW waveforms, we compute the re- 
sponse function of a space-based detector in heliocentric 
motion. Due to the motion of the LISA and eLISA/NGO 
constellations as they orbit (rotation and translation) it 
is more convenient to rewrite the response functions in 
terms of angles defined in a fixed Solar System Barycenter 
(SSB) coordinate system. The direction from the origin 
of the SSB reference frame to the origin of the EMRI 
reference frame is — n 



n = — (sin< 



sm v a sin c 



's- 



cos( 



i), (46) 



where (^s>^s) are spherical polar angles that determine 
the sky location of the EMRI with respect to the SSB 
frame. The relations between these angles and the angles 
(0(t),(j>(t)) that determine the sky location with respect 
to the detector reference frame are (see, e.g. [551 IM] 1: 

v^3 
2 

m 

where 



cos0(i) 



COS Wo — 



1 

2 

2nt/T- 



sin 9 S cos(2nt/T - </> s ) , (47) 
(48) 



= tan" 



•\/3 cos 9 S + sin 9 S cos (2-7rf /T — <j>$ ) 
2sin6> s sin (27rt/T - <£ s ) 



(49) 



and T = lyr is the period of the Earth's orbit around the 
Sun. The polarization angle ip describes the orientation 
of the "apparent ellipse" given by the projection of the 
orbit on the sky. It can be written in terms of (9§,(j) s ) 
and the angles describing the direction of the MBH spin 
with respect to the SSB reference frame (^K'^k)' 



tan-0 



jcos# K — V3 sin 9 K cos {2nt/T — </> K ) j 
-2 cos 9{t) {cos 9 K cos 9 S 



sm ( 



sm( 



i cos 



0} / 



sm tf K sm t/ s sm(0 K — <p s ) 

V3 cos (2irt/T) {cos 9 K sin 9 S sin <^ s 

- cos ^g sin 9 K sin cj) K } 

- \/3 sin (2irt /T) {cos #g sin 9 K cos </> K 

- cos# K sinflg cos^g} . (50) 
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In addition, the time of arrival of a gravitational wave- 
front at the SSB and at the detector will in general differ 
and are related by 

*ssb = *d + R sin6, s cos(27rf D /T - S ) - t% SB , (51) 

where R — 1 AU, i D is the time of arrival as seen in the 
detector reference frame, and tgg B is the initial time in 
the SSB reference frame: 



SSL! 



— It 



R sin S cos(2tt^ /T - 4> s ) . (52) 



This difference in arrival times gives rise to a Doppler 
modulation in the GW phase measured by LISA. To com- 
pute a waveform regularly sampled in time at the detec- 
tor, we need to generate a waveform unevenly sampled 
in the source frame (in which the time sampling is the 
same as at the SSB). This can be achieved employing the 
relations just introduced. 

The response of the detector to an incident GW can 
then be written as: 



Kit) 



V3 



F+{t)h + (t) + FZ(t)h x (t)} , (53) 



where a is an index for the different independent chan- 
nels of the detector. In the case of LISA we have two 
independent Michclson-likc interferometer channels that 
can be constructed from the LISA data stream and hence 
a = I, II. By contrast, cLISA/NGO will have only one 
independent channel (see Appendix [A] for a brief compar- 
ison of the detectors) and hence a = I for eLISA/NGO. 
The antenna pattern (response) functions, F^ ,x , are 
given by (see, e.g. [54]): 



F+ = - (1 + cos 2 9) cos(20) cos(2t/>) 
— cos 9 sin(20) sin(2^>) , 

Ff = - (1 + cos 2 9) cos(2(j>) cos(2t/>) 
+ cos 9 sin(20) sin(2^>) , 



F, 



i , - - (1 + cos 2 9) sin(20) cos(2^) 
+ cos 9 cos(20) sin(2^) , 



F£ = - (1 + cos 2 9) sin(20) sin(2V>) 
— cos 9 cos(20) cos(2?/>) . 



(54) 



(55) 



(56) 



(57) 



Here {9, <j>, ip) are as defined in Eqs. ( 47 1— ( 50 1 and spec- 
ify the sky location and orientation of the source in a 
detector-based coordinate system in terms of angles de- 
fined in a fixed SSB coordinate system. 



III. ELEMENTS OF SIGNAL ANALYSIS AND 
MODEL PARAMETER ESTIMATION 

The starting point for signal analysis is the detector 
data stream(s), s a . We assume that s a contains an 



EMRI GW signal, h a , and hence we can decompose it 
as 



»«(*) = M*) + n a {t), 



(58) 



where n a (t) is the noise in the detector, which we as- 
sume to be stationary, Gaussian and, in the case of LISA, 
that the two data streams are uncorrelated and the noise 
power spectral density is the same in each channel. Then, 
the Fourier components of the noise, which we denote 
with a tilde n a (f) ( see 1221 IM] for conventions on the 
Fourier transform that we use), satisfy 



1 



nM)nHf')) = ^ a0 5(f - f)S n (f) 



(59) 



where (•) denotes expectation value (ensemble average 
over all possible realizations of the noise), the asterisk 
now denotes complex conjugation, and S n (f) is the (one- 
sided) power spectral density (PSD) of the noise, which 
is given in Appendix [A] for both LISA and eLISA/NGO. 
The assumption of Gaussian noise means that the prob- 
ability of a particular realization of the noise n is given 
by 



p(n = n ) oc e 



-(n |n„)/2 



(60) 



where (-|-) denotes the natural inner product in the vector 
space of signals associated with the PSD S n (f) and is 
defined as 

W b)-^fdf '"^t^ , w 

a Jo b n {j) 

for any two signals a and b. The probability that a given 
GW signal h is present in a data stream s is thus 



p(s|h) oc e -(*-h|s-h)/2 



(62) 



The 'best-fit' waveform will be the one that maximizes 
(s|h) and, thus, it provides the maximum likelihood pa- 
rameter estimate. The expected signal-to-noise ratio 
(SNR), when filtering with the correct waveform, is 



SNR 



(h|h) 
rmsfhln) 



(63) 



where 'rms' stands for root mean square and the sec- 
ond inequality follows from the fact that the expectation 
value of (a|n)(b|n) is (a|b) [55]. In practice one considers 
a waveform template family that will depend on a set of 
parameters A, {h(t, A)}, and searches for the parameters 
that maximize the probability of a certain noise realiza- 
tion, i.e., the probability that a given waveform template 
is present in the data stream. Different realizations of the 
noise will lead to different values of the best-fit parame- 
ters. For large SNR, the best-fit parameters will follow 
a Gaussian distribution centered around the correct val- 
ues. Expanding exp(—(s — h|s — h)/2) around the best-fit 
parameters, A , by writing A = A + 5 A, we obtain the 
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following form for the probability distribution function 
for the errors 5\ 



lHt)\) = .\Vxp j -^T jk SX j SX k 



(64) 



where TV = v /dct(r/2?r) is the normalization factor and 



. %3 is the Fisher information matrix (FM) 



r jk - 



dh 



ah 

dX* 



(65) 



The variance-covariance matrix for the waveform param- 
eters is given by the inverse of the FM: 



{5X 3 5X k 



T- l ) jk [1 + 0(1/SNR)] , 



(66) 



and hence, we can estimate the precision with which we 
will be able to measure a particular parameter, A 1 , by 
computing the component T^ 1 of this inverse matrix, 
that is (see [57] for a detailed discussion): 

AX = J ((6X) 2 ) ~ Jr^ . (67) 



A. The Maximum-Mismatch Criterion 

Vallisneri |57j provided a consistency criterion to de- 
termine whether the SNR is high enough for the FM re- 
sults to be trustworthy, called the Maximum- Mismatch 
Criterion (MMC). The MMC criterion was suggested to 
assess when an estimation of the parameter errors based 
on a FM analysis would be reliable or not. Since the FM, 
, is built from the partial derivatives of the waveform 
template with respect to the parameters of the model, 
it can only represent the true GW signal, fe GW , cor- 
rectly if h(t, A) is linear in all the parameters, A, across 
a parameter space region of size comparable to the ex- 
pected parameter errors. This is the regime in which 
the Linearized- Signal Approximation (LSA) is valid. As 
we increase the SNR the errors become smaller and con- 
sequently the LSA is expected to work better. In the 
regime of validity of the LSA we can expand the wave- 
form template h(t, A) around the true source parameters, 
A tr , i.e. X 1 — X\ T + SX l with SX l being a small deviation 
in the parameters comparable with the parameter esti- 
mation error: 

h(t, A) = V + 5\* (dih) | Afr + ffih) \ xt + (68) 

Then, the likelihood [Eq. ( [62] )] can be approximated as: 
f (nln) ■ .(ahiah) 



+ 



SX j (a,h|n) } , 



(69) 



where the waveform template derivatives are evaluated 
at A = A tr . The applicability of the FM for parameter 



estimation is limited by the high-SNR requirement, in 
the sense that it can be a poor predictor of the amount 
of information obtained from waveforms depending on 
several parameters and detected with relatively low SNR. 
The MMC is given in terms of the ratio, r, of the LSA 
likelihood [Eq. (69)] to the exact likelihood [Eq. (60)]: 



2|logr| = [A\ i (d i h) xk - Dh AA J '(d,h 



'A* 



Dh (70) 



where Dh = h(A t fe r + AX k ) - h(A t fc r ) and AA* is the esti- 
mated error from the diagonal components of the inverse 
of the FM. The MMC is obtained by taking the maximum 
value of r over all parameters. 

The idea behind the MMC is to choose an iso- 
probability surface as predicted by the FM, and explore 
it to verify that the difference between the LSA and ex- 
act likelihoods is sufficiently small. Ratios, r, below some 
fiducial value are considered acceptable. If this condition 
is satisfied, we can believe that the FM is providing a 
reliable estimate of the parameter estimation errors. 



IV. PARAMETER ESTIMATION STUDIES: 
METHODS AND RESULTS 



In this section we describe the different techniques 
employed in our parameter estimation analysis and 
present the main results. We begin by charac- 
terizing the EMRI parameter space in our studies 
in DCSMG, {A*} 1=1 j Ar . There are 15 param- 
eters (the 14 of GR plus the DCSMG coupling 
parameter; see Table [IJ for a brief description): A = 
{M t , a, [a, CqjPq, iac Q, 0g, </jg, -Dl, "0q, Xq, '/'ol' 

where the subscript refers to the values of the 
corresponding quantities at the inspiral initial time. 

In order to simplify the computations involved in 
this study, we have restricted ourselves to a five- 
dimensional subset of the parameter space, given by 
A = {M„a, e X,D L /fi} (see Table [lj for their defini- 
tion). In this subset we have included those parameters 
that we have found to have the greatest correlation with 
the parameter £ = a ■ £ (see Table IT]), which controls 
the strength of the CS modifications (notice that in the 
MBH metric of Eqs. ([9])-(13l the CS parameter £ always 
appears multiplied by the spin parameter a and this has 
motivated the introduction of the combined parameter 
£) in a full parameter space investigation. We have also 
checked that the results we obtain do not change sig- 
nificantly when more parameters are added to the FM 
study. For the parameter estimation studies we consider 
two different EMRI systems, A and B, whose parameters 
are given in Table [TTJ These two types of systems differ 
in the values for the MBH mass, M % — 5 • 10 5 M Q for 
system A and M, — 10 6 M Q for system B. We fix the 
luminosity distance to D L = 1 Gpc, which roughly corre- 
sponds to the distance where we might expect the closest 
detectable sources to lie (see, e.g. [55] )■ Due to the fact 
that the inspiral time scales as ~ fi with the mass ratio, 
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TABLE I. Summary of the parameters that characterize an 
EMRI system in DCSMG. The angles (0 S ,0 S ) and (6> K ,0 K ) 
are spherical polar coordinates with respect to the ecliptic and 
the subindex stands for values of parameters computed at 
the initial time. The parameters with physical dimensions are 
indicated in square brackets. We set the luminosity distance 
to D h = 1 Gpc. 



Parameter 


Description 


M. 


MBH mass [M a ]. 


a = \S.\/M. 


MBH Spin [M.\. 


(i = m*/M. 


EMRI mass ratio. 


e 


Eccentricity of the particle orbit at to- 


Po 


Dimensionless semilatus rectum at t . 


e iac , 


Inclination of the orbit at t . 


C 


£-a [MS]. 


e s 


EMRI polar angle. 


4>s 


EMRI azimuthal angle. 


Ok 


MBH spin polar angle. 


<t>K 


MBH spin azimuthal angle. 


D L 


Distance from the SSB to the EMRI [Gpc]. 


V'o 


Angle variable for the radial motion. 


Xo 


Angle variable for the polar motion. 


4>o 


Boyer-Lindquist azimuthal angle. 



TABLE II. EMRI systems considered in the parameter esti- 
mation analysis. The table shows the values for the parame- 
ters that are considered in the FM computation (see Table [I] 
for the whole list of parameters). The rest of EMRI param- 
eters employed in our parameter estimation analysis are the 
same for both systems and their values are: m 4 = lOM©, 
inOj0 = 0.569, 6 S = S = 1.57, 9 K = K = 0.78, 4> = xo = 
0o = 0.78. 



System 


M. 


a/M. 


e 


C/M. 5 


D L / M [Gpc] 


A 


5- 10 5 


0.25 


0.25 


5- 10" 2 


5- 10 4 


B 


10 6 


0.25 


0.25 


5- 10" 2 


10 5 



the system A evolves faster than system B, which allows 
us to use smaller evolution times to obtain reliable results 
in that case. 

For these systems we evolve the trajectory using the 
geodesic equations given in Sec. |IIC| and the RR equa- 
Sec. ~ 

The ODEs that describe 



tions given in 
rithm outlined 



This is done using the algo- 
in Sec. IIIDJ 

geodesic motion are integrated for the angle variables 
[i{](t),x(t),<l>(t)] using the Bulirsch-Stoer extrapolation 
method [55] (see [501110] for details). The numerical code 
also contains routines which convert back and forth be- 
tween the different parameterizations of the orbit in DC- 
SMG; that compute the Cartesian orbital coordinates, 
velocities, and accelerations, the multiple moments; etc. 
The equations that evolve the constants of motion, E, L z 
and Q are integrated using simple finite difference rules. 
Then, we use the formulae of Sec. |II E| to compute the 
gravitational waveforms and the detector responses. 
In order to study how different the waveforms in DC- 



SMG are from GR, we have evolved our EMRI system 
during 0.5yr employing different values of the CS param- 
eter £ and the MBH spin a, and we have computed the 
following overlap function between a DCSMG and a GR 
waveform template: 



C [ h GR: h CSJ 



( h GR.I h Cs) 



GRl h GR.) ( h Csl h Cs) 



(71) 



0[h GR ,h GR ] = 



which is symmetric, O [h G 
and has the obvious property: 
O [h cs ,h cs ] = 1 . We also assume that the two wave- 
forms used for this overlap correspond to EMRIs with 
the same parameters, except for the CS parameter £ that 
vanishes for GR waveforms. The standard overlap de- 
fined in Eq. (611 has been computed using the FFTW 



library [61] for the Fourier transforms and simple integra- 
tion rules. We have computed this normalized overlap for 
a total of 121 EMRI systems which have 13 fixed paramc- 



= 0.25,p = 11, 



ters: M, = 5-10 5 M Q , m 4 = 10 M ( 
<2 mc = 0.569, S = 1.57, S = 1.57, 6 K = 0.329, 
cj> K = 0.78 , Djfi = 5 • 10 4 Gpc, % = 0.78 , Xo = 0.78 , 
4> a = 0.78, while the spin a/M t and the CS parameter 
£/Af 4 are varied in the interval [0,0.5] . The results are 
shown in Figure [2j where we can see how the projection 
of h cs onto h GR changes by modifying the values of the 
MBH spin a/M t and the CS parameter £/M 4 . In par- 
ticular, for higher values of a/M m and £/M 4 the overlap 
O [h GR , h cs ] decreases, since the difference in the evolu- 
tion of the SCO in GR and CS, produced by the dephas- 



ing introduced by the RR, increases [see Eqs. ( 16 )- ( 18 )] 
and, consequently, the deviations of h, 
enhanced. 



CS from h GR are 



CD 
> 

o 




FIG. 2. This 2D plot shows the symmetric normalized over- 
lap of Eq. ( 71 1 for EMRI systems with the following param- 



eters: M. = 5 • 10 5 M { 
<9 inc o = 0.569 



= 10 M , e = 0.25, po = 11, 



S = 1.57, 0s = 1.57 
5 ■ 10 4 Gpc, ^ = 0.78 , xo 



6 K = 0.329 , K = 0.78 , 
= 0.78, 0n = 0.78. The 



parameters f/M. and a/M. take values in the interval [0, 0.5] 
with a step of 0.05 for a total of 121 points. 
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We have also obtained the SNR in the frequency do- 



main using Eq. ( 63 ) . The computation of the FM requires 



the evaluation of the derivatives of the waveform tem- 
plates, c^h = dh./d\ l (actually of the response functions 
of the detector). Since the waveform templates/responses 
are generated numerically, the corresponding derivatives 
must also be evaluated numerically. For inner points in 
the EMRI parameter space (i.e. not near boundaries so 
that we do not need points outside the proper domains 
of definition of the parameters) we use the following five- 
point finite-difference rule: 

d,h = ^J— {h(A l + 26X 1 ) - h(X l - 28X 1 ) 



\28X l 
8 [h(A l 



<5A l ) - h(A l - 8X r )] } + O [(<5A l ) 4 ](?2) 

where S A* is the numerical offset in the parameter A 1 . For 
computations near the boundary or at the boundary of 
the parameter space we use instead non-centered finite- 
differences rules. Either the following three-point rule 

' {4h(A l + SX l ) - h(A ! + 2(5 A 4 ) - 3h(A 1 )} 



ah 



2 5X 
+ O [(5X) 2 } , 

or the following four-point rule 
1 



(73) 



ah: 



4 (5 A* 



{h{X l + 2SX*)-h(X t + 35X l ) 



+ 5 [h(A l + SX l ) - h(A 1 )] } + O [(6X 1 ) 3 ] . (74) 

It is known that computing numerical derivatives is a 
delicate task (see, e.g. [50 ). In the case of finite differ- 
ence formulas like Eq. ( |72[ ) the choice of the offset SX 1 is 
crucial. An offset too small will produce high-order can- 
cellations in the numerator beyond machine precision. In 
contrast, an offset too big may mean higher order terms 
in the Taylor series expansion of the waveform become 
important. In both cases we will be far from a reasonable 
approximation. Therefore, we have done investigations 
that survey wide ranges for SX l in order to find intervals 
where the derivatives have good convergence properties. 

Once we have obtained a FM, r^, that converges in 
a certain range of offsets SX 1 , we estimate the expected 
measurement error in the parameters by using Eq. ( 67 1 . 



Since FMs for EMRI waveforms have very large condition 
numbers (the ratio of the largest to the smallest eigenval- 
ues), we use an LU decomposition to invert them, writing 
the matrix as the product of a lower triangular matrix 
and an upper triangular matrix In addition, to as- 
sess whether the error estimates obtained arc reliable or 



not, we use the MMC defined in Eq. (70). We have eval 



uated the MMC criterion for all the results presented in 
this paper and, unless otherwise specified, they satisfy 
this criterion with values of | logr| ranging from 10~ 4 to 
0.5. 

We have stated before that RR effects change the rel- 
ative phase between waveforms in DCSMG with respect 
to GR. Now, we are going to study their impact on pa- 
rameter estimation. Firstly, we will compare the param- 
eter estimation errors for systems evolved under RR and 



systems that do not radiate, preserving the constants of 
motion, i.e. always have the same orbital parameters. 
These results have been obtained assuming that the de- 
tector is LISA. At the end of this section we present 
some results for eLISA/NGO. The results for system A 
with and without RR and for different evolution times, 
0.3 , 0.5 , 1 yr, are shown in Table IIIII The 



T. 



cvol 



0.1. 



upper part of the table contains the results for the evo- 
lutions with RR and the lower part of the table shows 
the results without RR. In both cases we also show the 
value of the MMC test, i.e. the quantity | log r | defined 
in Eq. (70). We do not show results for T evol = 0.1 yr 



without RR in Table IIIII since we did not obtain reli- 
able results (according to the MMC criterion). From 
these results and others we have obtained for other sim- 
ilar EMRI systems we can say that the typical measure- 
ment accuracies for the five most important parameters 



" 6 M. 

■ 10" 



io- 



are: AlogM. ~ 10~ 3 , Aa ~ 10" 

A log C ~ 10~ 2 and Alog(D L //x) ~ 10~ 2 . Comparing 
the results in Table III we see that the inclusion of the 
RR improves the SNR of the signals. It also improves 
the parameter estimates, in particular those of the spin, 
a, and of the CS parameter, £. This is partially due to 
the increase of the overall SNR due to RR, but even after 
rescaling to a fixed reference SNR we see an improvement 
in the parameter measurement accuracies when RR is in- 
cluded. As one could expect due to the adiabatic nature 
of the RR (see e.g [S3]), the improvement with the inclu- 
sion of RR is more significant for longer evolution times. 

We can also explore how the error estimates change 
with the spin parameter a. To that end, we did simula- 
tions of systems A and B using the following values of the 
spin: a/M. = 0.1 (A^Bj), a/M. = 0.25 (A 2 ,B 2 ), and 
a/M. = 0.5 (A 3 ,B 3 ). The initial semilatus rectum was 
set to p Q = 11M., which means that after evolving sys- 
tems A l -A 3 for a total time of T evol = 0.5 yr and systems 
B 1 -B 3 for a total time T ovol = 1.5 yr, the final semila- 
tus rectum, pp is approximately 8 M. for all systems. 
The parameter estimation errors are shown in Table IV 
(the upper part corresponds to simulations of system A 
whereas the lower part corresponds to simulations of sys- 
tem B). The first thing that we notice is that the smaller 
the spin parameter a/M. becomes, the better the pa- 
rameter estimate for the CS parameter Q. In particular 
AC ~ 2.8 • lO- 2 for system A l and AC ~ 1.4 • 10" 2 for 
system B 1 . The reason for this is quite simple. The CS 
modifications affect a single MBH metric component, g t ^ 



[Eq. (13)], which contains the CS parameter £/Mf, mul- 
tiplied by the spin parameter a. The unperturbed Kerr 
metric component is proportional to a, and so the rela- 
tive change in this metric coefficient due to the addition 
of the DCSMG correction is proportional to £. Since we 
keep C = a £, fixed as we vary a, the value of £ increases as 
a decreases and so the CS correction to the MBH metric 
is larger relative to the leading order Kerr metric term. 

The values of the SNR that we obtain for systems A Y , 
A 2 , and A 3 are, 48.1, 46.5, and 44.7 respectively. In the 
case of systems B 1: B 2 , and B 3 , the values of the SNR are 
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TABLE III. Error estimates for LISA and the EMRI system A (see Table [Tip using RR (upper part of the table) and without 
using RR (lower part of the table). Each column contains the estimations for a given evolution time (T cvol = 0.1 , 0.3 , 0.5 , 1 
yr) and shows the corresponding SNR of the EMRI signal. 



With RR 


T — 

evol 


0.1 yr 




T — 

evol 


0.3 yr 




T — 

evol 


0.5 yr 




T — 

evol 


lyr 








SNR= 


14.5 




SNR= 


43.2 




SNR= 


55.4 




SNR= 


73.5 






X 1 


AA ! 


1 logr| 




AA l 


1 logr| 




AA l 


|logr| 




AA l 


lo 


gr| 




logM. 


1.4- 10" 1 


1.5- 10" 


l 


9.2 • 10 -3 


1.1 • 10" 


i 


4.5 ■ 10~ 3 


1.5- 10" 


-l 


9.3- 10" 4 


2.4- 


10" 


1 


a/M. 


1.2 • 10~ 4 


3.4- 10 - 


i 


1.5- 10" 5 


2.0- 10- 


i 


4.9 • 10~ 6 


5.3- 10- 


-2 


1.5 • 10" 6 


2.3 • 


10" 


1 


e 


5.2 • 10" 6 


5.2- 10" 


2 


9.6- lO" 7 


3.0- 10" 


2 


5.0- 10" 7 


9.7- 10" 


-3 


2.8 ■ 10 -7 


6.0- 


10- 


3 


logC 


1.1 


9.3- 10" 


1 


1.5 • 10" 1 


3.1 • 10" 


1 


4.9 • 10~ 2 


3.1 ■ 10- 


-2 


2.0 ■ 10 -2 


1.5 • 


10- 


1 


log(A» 


2.0- 10" 1 


2.0- 10- 


1 


1.5- lO" 1 


4.1 • 10" 


4 


1.8 ■ 10 -2 


1.6- 10- 


-4 


1.3 ■ 10 -2 


2.6- 


10" 


4 


With no RR 








T — 

evol 


0.3 yr 




T — 

evol 


0.5 yr 




T — 

evol 


1 yr 














SNR= 


38.4 




SNR= 


46.8 




SNR= 


54.6 







y 


AX' 


logr| 


AX' 


logr| 


AA* 


logr| 


logM. 


8.3- lO" 3 


4.3 • 10~ 2 


3.9- 10" 3 


3.7- lO" 2 


6.6- 10" 4 


2.3- lO" 4 


a/M. 


2.3- lO" 5 


3.3- 10" 1 


1.4- 10 -5 


2.4 ■ 10" 1 


7.4- 10 -6 


1.6- 10" 1 


e 


1.0- 10" 6 


5.0- lO" 2 


6.7- 10~ 7 


3.6- lO" 2 


1.4 • 10~ 6 


1.6- lO" 3 


logC 


2.5- lO" 1 


6.1 • 10" 1 


1.5- lO" 1 


3.6- 10" 1 


1.1 • 10" 1 


3.2 • lO" 1 


log (Djfj.) 


2.8- lO" 2 


4.8 • 1Q- 4 


2.1 • lO" 2 


4.8 • lO" 4 


1.8 • 10" 2 


3.8- 1Q- 4 



50, 46, and 49. Notice that the SNR varies by modifying 
the value of the spin parameter a, albeit not very much. 
This dependence of the SNR on the system parameters 
is expected in the region of the parameter space where 
the FM can be linearized (see e.g. [64] and [65]). 

Overall, the parameter estimation errors for system A 
have the magnitudes: 

A log M. ~ 5 • 1(T 3 , Aa ~ 5 • 1CT 6 M. , (75) 
Ae - 3 • 1CT 7 , A log C - 4 • 1CT 2 , (76) 

Alog(^ L / A/ )-2-10- 2 . (77) 

In the case of system B they are: 

A log M. - 6 • 1(T 4 , Aa - 3 • 1(T 6 M. , (78) 
Ae - 1CT 7 , A log C ~ 2 • 1CT 2 , (79) 

Alog(i? L / Ai )^2-10- 2 . (80) 

The order of magnitude is roughly the same for both 
systems, but in general the estimations for system B are 
better than those for system A, since the MBH mass for 
system B is larger than the one for system A and the 
integration time is longer, so there are more observed 
waveform cycles. 

The parameter error estimates presented are for a fixed 
value of the parameter £ = £-a. Since the spin parameter 
a/M, is fixed and is the same for both systems A and B 
in Table [Tl] this means that in the previous results the 
CS parameter £/AT 4 was fixed. Now, we present results 
for the EMRI system A for different values of the CS 
parameter £/Mf. We have considered the following par- 
ticular values: f = 0.05M. 4 , £ = 0.1M. 4 and £ = 0.2M. 4 . 
The results obtained for the estimation of the parameter 
errors of A = {M # , a/M,, e ,p , D L //j,} are shown in 
Table [Vj Due to the fact that the dependence on £/M 4 



and on a/M. are different in the MBH metric compo- 
nents and in the evolution equations (see Sec. [TT| , one 
would expect a different dependence of the error esti- 
mates when varying both parameters independently. By 
comparing the results of Tables |IV| and [V] we can see 
that modifying the value of the CS parameter £ /M 4 only 
affects significantly the error estimate of the CS parame- 
ter itself, whereas modifying a / M, has a significant effect 
on the error estimates of all the parameters employed in 
our study, and in particular on £. 

Up to now, all the parameter estimation results pre- 
sented refer to the LISA detector. We now present some 
results for eLISA/NGO [6 . In order to more easily 
compare with the results obtained for LISA, we nor- 
malize to a fixed SNR, since the SNR for eLISA/NGO 
is around two times smaller. We considered system A 
with three different values of the spin parameter a/M. 
namely a/M, = 0.1, 0.25, and 0.5. The results ob- 
tained are quoted in Ta ble |VI| Comparing them with 
the ones quoted in Table |IV| for LISA, we can see that 
the parameter estimation accuracy does not change ap- 
preciably when the noise curve of LISA is changed for 
the one of eLISA/NGO and so, all previous results can 
be considered to apply to eLISA/NGO as well, with the 
corresponding SNR corrections. 



V. PLACING A BOUND ON THE CS 
PARAMETER 

One application of the framework we have developed to 
perform parameter error studies in DCSMG is to try to 
put bounds on the CS parameter £, which is the combina- 
tion of CS coupling constants and the gravitational con- 
stant that controls deviations from GR in the dynamics 
of EMRIs. This question has already been investigated in 
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TABLE IV. Error estimates for LISA and the EMRI systems A and B. The results shown have been obtained for different 
values of the initial eccentricity, eo, and MBH spin, a. The evolution time for these systems is: T ovol = 0.5 yr (system A) and 
T ovol = 1.5yr (system B). The superscript "'" on a given result indicates that the corresponding Fisher Matrix did not satisfy 
the MMC criterion, nevertheless we include the results for the sake of completeness. 



system A 


e = 0.1 


e = 0.25 


e = 0.5 


a/M. 




0.1 




0.25 


0.5 T 






0.1 




0.25 




0.5 






0.1 




0.25 




0.5 




logM. 


4.2 


• 1(T 


a 


4.1 • 10 _a 


3.0 ■ 10" 


a 


3.7 


• io- 


3 


4.3 • 10-' J 


4.4 


■ io- 


a 


1.4 


■ io- 


a 


5.0- 1Q- J 


4.9 


• io- 


3 


a/M. 


5.0 


• io- 


5 


6.0- IO" 6 


8.0- 10- 


6 


3.2 


■ 10" 


6 


5.2 • IO -6 


7.2 


• 10" 


6 


4.0 


• io- 


6 


4.6 • IO -6 


6.0 


■ io- 


6 


e 


2.3 


• io- 


(i 


2.4 ■ IO' 6 


1.4- 10" 


6 


4.9 


■ io- 


-7 


8.6- IO" 7 


9.2 


■ io- 


7 


2.0 


■ io- 


7 


3.3- 10" 7 


3.3 


• io- 


-7 


logC 


7.5 


• 10" 


2 


9.9- 10" 2 


9.0- 10" 


2 


2.8 


■ 10" 


2 


4.9 • 10~ 2 


6.6 


• 10" 


2 


5.1 


• 10" 


2 


3.5- IO" 2 


4.3 


■ 10" 


2 


log (DJu) 


1.9 


• io- 


2 


2.0 ■ 10 -2 


2.1 ■ 10" 


2 


1.7 


• 10- 


-2 


2.0 • 10 -2 


2.1 


■ 10- 


2 


1.9 


■ 10- 


2 


2.3 • IO -2 


2.4 


• 10- 


-2 


system B 


e = 0.1 


e = 0.25 


e = 0.5 


a/M. 




0.1 




0.25 


0.5 T 






0.1 




0.25 




0.5 






0.1 




0.25 




0.5 




logAf. 


1.1 


• 10" 


a 


1.1 ■ 10 _a 


5.3- 10" 


4 


8.7 


■ io- 


4 


9.0 -10" 4 


9.9 


• 10" 


4 


6.1 


■ io- 


4 


6.4- 10 -4 


6.8 


• io- 


4 


a/M. 


4.8 


• io- 


(i 


7.7- IO" 6 


5.1 • 10" 


6 


3.2 


■ io- 


-6 


4.3 • 10~ 6 


5.2 


• io- 


6 


1.7 


• io- 


6 


2.1 • IO" 6 


2.7 


■ io- 


-6 


e 


1.5 


• 10" 


6 


1.9 ■ IO -6 


6.0- 10" 


7 


4.1 


■ 10" 


7 


4.3 • IO -7 


4.5 


■ io- 


7 


9.6 


■ 10" 


8 


1.0- 10" 7 


1.1 


• 10" 


7 


logC 


6.4 


• HT 


2 


8.5- 10~ 2 


7.0- 10" 


2 


3.7 


■ io- 


-2 


4.3 • IO -2 


5.3 


• 10" 


2 


1.3 


• io- 


2 


1.6 • IO -2 


2.1 


■ io- 


-2 


log {DJ») 


2.6 


• 10" 


2 


2.7- 10~ 2 


2.0- 10" 


2 


2.3 


• io- 


2 


2.4- IO -2 


2.6 


• 10- 


2 


2.0 


■ 10- 


2 


2.1 • 10~ 2 


2.2 


• 10- 


2 



TABLE V. Error estimates for LISA and the EMRI system 
A in Table [II] obtained by changing the value of the CS pa- 
rameter £. As we can see, by increasing the value of the CS 
parameter, f, its error estimate, Alog£, improves, whereas 
the rest of the error estimates remain roughly constant. 





i/Mt 


= 0.05 


C/A/. 4 = 


0.1 


i/Mt = 0.2 


logAf. 


4.4- 




4.2 • 10 


-3 


4.5 ■ 10 _a 


a/M. 


4.9- 


10" 6 


4.7- 10 


-6 


4.9 • 10~ 6 


e 


4.9 • 


io- 7 


4.9 ■ 10 


-7 


5.0- IO" 7 


logC 


1.9- 


io- 1 


9.5- 10 


-2 


4.9 • 10" 2 


log (D L /n) 


1.8- 


10 -2 


1.8- 10 


-2 


1.8 ■ 10 -2 



TABLE VI. Parameter estimation results for eLISA/NGO 
and the EMRI System A. evolved for T cvol = 2 yr. The cor- 
responding SNR is ~ 15 . The value of the initial eccentricity 



is e n = 0.25 . 





a/M. = 0.1 


a/M. 


= 0.25 


a/M. = 0.5 


logAf. 


9.0 • 10" 4 


9.8- 


10 -4 


1.3- 10 _:i 


a/M. 


3.2 • 10" 6 


2.8- 


IO" 6 


3.9- 10" 6 


e 


5.1 ■ 10" 7 


5.2- 


io- 7 


5.7- 10~ 7 


logC 


6.0- IO" 2 


7.3- 


IO" 2 


9.6- 10~ 2 


log(f? £ /At) 


6.4- 10" 2 


7.0- 


IO -2 


7.5- IO" 2 



the literature, but using astrophysical systems different 
from EMRIs. 

In the case of non-dynamical CS gravity, although 
the scalar field i9 has no evolution equation, it can 
be prescribed a certain time evolution that has an as- 
sociated time-derivative § and timescale, r cg = l/i?. 
Bounds are normally written in terms of constraints 
on £ 2 /t cs , where H 2 is the characteristic length scale 
and equals the coupling constant a introduced earlier. 
Strong bounds on this combination were first obtained 
by Yunes and Spergel [66) but refinements introduced 



by Ali-Haimoud [ST] set the bound to 0.2 km, which is 
three orders of magnitude better than the Solar System 
bound [55], obtained from data from the LAGEOS satel- 
lites orbiting the Earth [69 . 

For dynamical Chern-Simons gravity, which we con- 
sider here, the bound is normally expressed as a bound 
on C 1 / 4 . The first bound was quoted by Yunes and Preto- 
rius |25j and was < 10 4 km. However, in a recent pa- 
per by Ali-Haimoud and Chen [35] they took into account 
the fact that the CS solution for the spacetime outside a 
rotating star is not the same as that outside a rotating 
black hole, and also that the CS correction can only lead 
to a decrease in frame-dragging effects and thus cannot 
be constrained by an upper bound on the precession, but 
only by a positive lower bound that lies below the GR 
value. The CS-induced precession that gives the bound 
quoted in Yunes and Pretorius is two orders of magnitude 
larger than the GR precession, which means that bound 
can probably not be trusted. Ali-Haimoud and Chen [35] 
argue that the best constraints at present are therefore 
those that come from Solar System measurements, based 
on data from the Gravity Probe B satellite [70] and also 
from the LAGEOS satellites [69 , which are much weaker. 
The bound that they get is then £ x / 4 < 10 8 km. In this 
paper we compare our results with this weaker but more 
robust bound. 

The basis for the computation of our bound is the fol- 
lowing. We assume that GR is the correct theory to 
describe EMRI dynamics and hence assume that mea- 
surements made by LISA are compatible with £ = 0. 
Then, by estimating the error on the measurements of 
£, A£ [obtained using Eq. (67)], we can set a bound of 
the following type: £ < A£ . Different EMRI systems 
will provide different constraints. But since £ is a uni- 
versal quantity, in particular the same for all EMRIs, 
we just need to look for the EMRI system that provides 
the best constraint. We have performed several compu- 
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TABLE VII. Parameter estimation errors for EMRI systems in GR (i.e. on the A-parameter surface determined by £ = ( — 0). 
The parameters common to all these systems are: M. = 5 • 10 5 M Q , /i = 2 • 10~ 5 , 9o, ~ 4>s = 1.57 rad, K = 0.392 rad, 
4> K = 0.78 rad, D h /fi = 5 ■ 10 4 Gpc, and ipo = Xo = 4>o — 0.78 rad. The last two columns contain the error estimate for the CS 
parameter £ and the MMC criterium figure of merit associated with the CS parameter £ . 



a/M. 


e 


Po 


e ioo , 


T cvo i (yrs) 


A£/M 4 


1 log r | 


0.5 


0.7 


10.0 


0.15 


0.5 


5.76 • 10~ 8 


0.91 


0.45 


0.7 


10.0 


0.15 


0.5 


1.86- 10~ 7 


0.84 


0.5 


0.7 


10.0 


0.15 


0.5 


6.23 • 10 -8 


0.89 


0.5 


0.85 


11.0 


0.15 


0.5 


6.20 • 10" 8 


0.7 


0.5 


0.85 


11.0 


0.15 


1.0 


6.10 • 10" 8 


0.57 



tations with EMRI systems whose common parameters 



are: M = 5 • 10 5 M 



o 



2 • 1(T 5 , 



1.57 rad, 



9 K = 0.392 rad, cj) K = 0.78 rad, DJu = 5 • 10 4 Gpc, and 
— Xo _ — 0.78 rad. We show some relevant results 
in Table VII for EMRIs with spin parameter a/M. ps 0.5 



(the rest of parameters can be found in the caption of 
the table). Since we are differentiating about zero, the 
numerical evaluation of the £ derivatives must be per- 
formed using a one-sided derivative. In particular, we 
have double-checked some of these results using both the 
3-point and 4-point rules given by Eqs. (73 1 and (74) 
respectively. We note that, even though we are using 
system A (M. = 5 • 1O 5 M ) for this study, the values ob- 
tained for the MMC with T ovo i = 0.5 were slightly above 
the reference threshold of 0.5 that we used elsewhere in 
our study. This fact could be connected to using the 
one-sided derivative in our calculations. 

From the error estimates for the £ parameter shown in 
Table IVlTl, we find 



A£/M 4 < 10~ 7 . 
which, in suitable units, becomes 

£ 1/4 < 1.4- 10 4 km. 



11) 



(82) 



This result, a prediction for LISA measurements, is al- 
most four orders of magnitude better than the bound 
£i/4 < 10 8 km 

given in [32 . The corresponding esti- 
mate for eLISA/NGO can be found be rescaling the £ 
constraint by the SNR, but since the bound scales only 
as the one-fourth power, the bound for eLISA/NGO is 
essentially the same. 



VI. CONCLUSIONS AND DISCUSSION 

In this paper we have examined how well a space-based 
GW detector like LISA or eLISA/NGO can discriminate 
between an EMRI system in GR and one occurring in a 
modified gravity theory like DCSMG. To do this, we have 
extended previous work in I27j by introducing two key 
components. The first is the inclusion of RR effects driv- 
ing the inspiral. We have constructed a waveform tem- 
plate model using an adiabatic-radiative approximation 



following the NK waveform model [31] . In this approxi- 
mation, the inspiral trajectory is modeled as a sequence 
of geodesies whose constants of motion are updated using 
formulae for the fluxes of energy, ^-component of the an- 
gular momentum, and Carter constant, that were derived 
for general relativistic inspirals in |29j using a combina- 
tion of PN approximations and fits to results from the 
Teukolsky formalism. 

The second key improvement made in this paper 
is the use of the Fisher matrix formalism to esti- 
mate errors in parameter measurements. We have 
explored a five-dimensional subspace of the fifteen- 
dimensional parameter space of EMRIs in DCSMG (see 
Table IB. The parameters that span this subspace are 
{M.,a/M 9 ,e ,£,D L /fj,}. We have focused our studies 
on two types of systems (see Table [TTJ) , with masses 
1OM + 5 ■ 10 5 M Q and 10M Q + 1O 6 M . The parameter 
error estimates are summarized in Eqs. (77 1 and (80 1. 



For both systems, and assuming a LISA detector, we es- 
timated the measurement error on the logarithm of the 
CS parameter £ as Alog£ ~ 10~ 2 . Therefore, a space- 
based GW detector like LISA should be able to discrim- 
inate between GR and DCSMG. In the case that DC- 
SMG is the correct theory describing the strong gravita- 
tional regimes involved in EMRI dynamics, such a detec- 
tor should be able to provide a good estimation of the 
CS parameter that controls the deviations from GR. We 
have also explored how these parameter error estimates 
change with the spin parameter a/M. (see Table IV I. We 
have found that by decreasing a/M. the parameter mea- 
surement precision of the CS parameter improves, while 
modifying the value of £ (see Table [Vj does not affect 
significantly the precision of parameter estimates for the 
range of other system parameters included in our analy- 
sis. 

For the case of eLISA/NGO we have presented some 
parameter error estimation results for system A in Ta- 
ble |nj In order to compare with LISA results we have 
normalized these results to a fixed SNR (the eLISA/NGO 
SNR is approximately a factor of two smaller than the 
LISA one) . Results for three values of the spin parameter 
(a/M. = 0.1 , 0.25 , and 0.5 ) are given in Table [Vl] The 
conclusion is that the parameter estimation accuracy at 
fixed SNR does not change significantly relative to the 
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LISA results. The LISA results can therefore be applied 
to eLISA/NGO by applying the appropriate SNR correc- 
tion. 

Finally, we have used our parameter estimation frame- 
work to put bounds on the CS parameter £. By assuming 
that GR is the correct theory of gravity we have found 
that LISA could place a bound £ 1/4 < 1.4 • 10 4 km, which 
is almost four orders of magnitude better than the bound 
obtained in [35] using Solar System data. 

The results presented in this paper can be extended in 
a number of ways by adding more elements to the wave- 
form model that we employ. For example it can be done 
by: (i) Using a higher order approximation for the MBH 
geometry in DCSMG; (ii) including CS corrections to the 
RR formulae, in particular to introduce the effects of the 
CS scalar field in the RR mechanism; (iii) adding more 
multipole moments to the gravitational wave expansion 
formulae; etc. In addition, we have focussed our study on 
a few EMRI systems, so it would be useful to carry out a 
more exhaustive study of the parameter space, although 
this would be a costly task in terms of computational 
resources. Such extensions to the present work would 
allow us to consider systems that might be of greater in- 
terest from the point of view of improving the parameter 
estimation results. The approximations underlying our 
model prevent us from considering systems with spins 
higher than ajM m = 0.5 and strong CS couplings. How- 
ever, a better search of the parameter space would allow 
us to identify systems that provide the best parameter 
estimates and the strongest bounds on the CS parameter 
t 

There are other extensions of this work that are also 
interesting. In particular, it would be useful to assess the 
systematic errors that would arise if GR waveform tem- 
plates were used to detect EMRIs that are actually de- 
scribed by DCSMG. This could be done using the formal- 
ism developed by Cutler and Vallisneri [7T] to estimate 
systematic errors that arise from model uncertainties. Fi- 
nally, we could apply some of the tools and techniques 
used in the present work to study other modifications of 
gravity, different from the CS correction and in this way 
to exploit the potential of the connection between gravi- 
tational wave astronomy and high-energy physics [721. 
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Appendix A: LISA and eLISA Power Spectral 
Densities 

In this paper we assume that the GW detector is either 
LISA @1 [5J or eLISA/NGO 0. LISA is a space-based 
GW detector concept that consists of a quasi-equilateral 
triangular constellation of three identical spacecrafts 
with an inter-spacecraft distance of L = 5 ■ 10 9 m. Each 
spacecraft follows a heliocentric orbit that trails behind 
the Earth at a distance of 5 • 10 10 m (equivalent to 20 de- 
grees) in such a way that the LISA constellation faces the 
Sun, slanting at 60 degrees to the ecliptic plane. These 
particular heliocentric orbits were chosen such that the 
triangular formation is maintained throughout the year, 
with the triangle appearing to rotate about the center of 
the formation once per year. Each spacecraft contains 
two free-falling test masses whose distance is monitored 
by 6 laser links. In contrast, the eLISA/NGO constella- 
tion has an inter-spacecraft distance of L = 10 9 m. More- 
over, only one of the spacecrafts will contain two free 
falling masses and service two arms of the constellation, 
while the other two will have only one proof mass and 
service one arm. This effectively reduces the detector re- 
sponse from having two independent Michelson channels 
to just one. 

An essential ingredient required in the detector re- 
sponse is a model for the noise affecting the obser- 
vations. This may be described in terms of the 
one-sided noise power spectral density, S n (f). For 
LISA, this has three contributions: instrumental noise, 
S„ st (f), confusion noise from short-period galactic bi- 
naries, <S^ al (/), and confusion noise from extragalactic 
binaries, S£ xgal (/) [51]: 

S n = mm{ST + S^ al , Sr* + 1 + S^}(M) 

where the different noise contributions are given by: 

S™\f) = exp [nTj ssion ^ (9.18 x W^f^ 

+ 1.59 x 10~ 41 + 9.18 x l(r 3S f 2 ) Hz" 1 , (A2) 

/ f \-7/3 

ST\f) = 2.1 x 10- 45 (^J Hz" 1 , (A3) 

S^\f) = 4.2 x 10- 47 ( J^) V3 Hz- 1 , (A4) 

with dN/df the number density of galactic white dwarf 
binaries per unit frequency, T^ssion the lifetime of the 
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LISA mission, and k the average number of frequency 
bins that are lost when each galactic binary is fitted out. 
The particular values that we use correspond to (see, 
e.g. [51] ) K ps 4.5 and: 



d7 = 2x10 -(t-J 



11/3 



(A5) 



The noise curve for eLISA/NGO is given by [S]: 



S n (f) = 



AQ _d_ Q _|_ C 

acc *^gn t u omn 

L 2 



0.205c 



(A6) 



where 5„ 



5 sn and .!5 omn 



are, respectively, the power 



spectral density of the residual acceleration of the test 
masses, of the shot noise and of other measurement 



noises. These are given by: 





= 1.37 


10" 


32 L + 


s s M) 


= 5.25 


10" 


23 m 2 Hz 


q 

omn 


= 6.28 


10" 


23 m 2 Hz 



10- 4 Hz\ Hz 2tt _ 1/a . 

^m 2 Hz XA7) 



f 4 



(A8) 
(A9) 



Appendix B: Evolution of the Constants of Motion 

In [2H1I3I]j to compute EMRIs in general relativity, the 
fluxes on the right-hand sides of Eqs. ( 31 )-( 33 ) were spec- 
ified by approximate, weak-field, formulae, augmented 
with corrections to ensure the behavior was not patholog- 
ical for near-circular or near-polar orbits and augmented 
by fits to numerical solutions of the Teukolsky equation. 
These formulae look as follows (note that in this paper 
we are using a dimensionless semilatus rectum instead of 
the semilatus rectum of which has units of A/.): 



dE 
~dt 



dt 

dQ 
dt 



(1 - e 2 ) 3 / 2 
N 5 (p,L 



(1 _ e 2 )-^ 2 (E) 2PN (p, i, e, a) - (E) 2PN (p, t, 0, a) 



N 4 (p, t) 



(Q^pnG^.M) 



(1 - e 2 f' 2 (1 - e 2 r 3/2 (4) 2 PN(P, h e, a) - (L z ) 2PN (p, t, 0, a) + (L,) fit 



= (l-e 2 ) 3/ V%^,«) 



(l_e 2 )" 3 / 2 



Q 



Q 



(p,L,e,a) - 



Q 



2PN 



C7 



(p,t,0,o) 



2PN 



2 taut ^ (Lj fit + — 2 (t) fit 

sin t 



(Bl) 
(B2) 



(B3) 



where the coefficients N/s are 



N 1 { P ,l) = P M: 



N 4 ( P ,L) = pMl 



M 2 



pE{p 2 +q 2 )-2q[^--qE 



L. 



(2-P) 



2 9 £; 



(B4) 

3irc 

(B5) 
(B6) 



where the subscript 'circ' indicated that these coefficients 
are evaluated for a circular orbit defined by the argu- 
ments p and l. In these expressions, q denotes the di- 



mensionless spin parameter of the MBH 
a 



M. 



< q < 1 



(B7) 



In Eqs. (B1|-(B3), the fluxes (E) 



^2PN' 



z)2PN> 



and 

(Q) 2PN are the 2PN approximations to the averaged evo- 
lution of the energy, angular momentum in the spin direc- 
tion, and Carter constant. They are modifications of the 
original expressions given in |73j but corrected to avoid 
unphysical features that they exhibit for nearly circular 
(e ps 0) and for nearly polar (t ps tt/2) inspirals. The 
corrected 2PN fluxes have the following form 
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(^)2PN — 



32 ml (1 - e 2 ) 3 / 2 



(4) 



2PN 



(Q) 



2PN 



5 Ml p 5 




'527 6533 2 


\ q 2 


v 96 + 192 6 


) p 2 


32 m 2 (1 - e 2 ) 


3/2 


5 M. p 7 / 2 




|^ 5l2 (e)cost 


1 


64 m 2 (1 - e 2 ) 


3/2 



3i( e ) - pYj2 -92 ( e ) cos ^ - - 9 3 (e) + s 4 (e) - ^ g 5 (e) + ^ <? 6 (e) 



sin 2 t 



flr 9 (e) cosi + -|- |g?o( e ) - .9w( e ) cos 2 /.} - -.9 n (e) cos/, 



45 37 



. 2 

sm t 



5 M. 



,7/2 
2 
2 



4 5i 3 (e) + Is ( fi4( e ) ~ sin2 1 



Q sin i 
45 



fg( e ) - -4/2 5 io( e ) cost ^ £.9n( e ) + ^72 5i 2 ( e ) 



1 

P 



7T 

g372 



(B8) 



(B9) 



(BIO) 



I 

where the various e-dependent coefficients are 



73 2 37 4 i . 73 823 2 949 4 491 fi . . 1247 9181 2 

o 1 (e) = H e H e 4 , g,(e) = 1 e H e H e° , g,(e) = 1 e , 

yiV y 24 96 ' y2W 12 24 32 192 ' y3V ; 336 672 ' 

. , , 1375 2 , , 44711 172157 , . . 33 359 2 

■ 94(e) = 4+ ~W e ' ffs(e) = M72" + ^592" e ' 9e{6) = 16 + ~S2 e ' 

, . 8191 44531 2 , , 3749 5143 2 , s , 7 2 

9 7 e = 1 e , Qo ( e ) = e , g Q ( e ) = 1 H — e , 

J7V ; 672 336 ' ysK ' 336 168 ' y9W 8 

a , , 61 63 2 95 4 h . , 61 91 2 461 . , , 1247 425 2 

5l0 (e) = ^4 + y e +64 ' 5lo(6) = T + T 6 + "64" 6 ' 5ll(e) = 336 + 336 6 ' 

/ x „ 97 2 , , 44711 302893 2 , . 33 95 2 

o 1 _(e)=4H e , q.Je) = 1 e 2 , 914(e) = 1 e . (Bll) 

ym y t- g , yi 3 v ; gQ72 t 6Q48 , y 14 y ' m ^ w V > 



The equation for the evolution for the Carter constant 
has an additional improvement with respect to the one 
in |73j . where a simple but accurate prescription for the 
Carter constant was given by assuming that the incli- 
nation angle evolution due to GW emission is negligible 
(see HI] for supporting evidence of this). That is, 



I fa leads to Q w 2(LJL Z )Q via Eq. (j30j). The im- 
provement introduced in |29j consists of adding the next- 
order spin-dependent PN correction. 

The final ingredient comes by adding fitting func- 
tions to the results of Teukolsky-based computations 
for circular- inclined orbits [30,. The expressions for 
the Teukolsky fitted fluxes (to data provided by Scott 
Hughes) are 
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Similarly, a good fit to the evolution of i is given by 
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where the values of the numerical fitting coefficients are: 



d\ = -10.7420, d\ 
d c 2 = -33.7090 , 



* _ 



28.5942 , 
c\ = -28.1517, 



4 2 = 2.37258, 
4 = -7.61034, 
4 = 306.119, 
4 7 = 224.227, = - 

c% = -0.645000 , 
/ 2 ° = 483.266 , 
fl = 82.0780 , 



4 = -66.6584, 
c\ = 128.878 , 
4 = 40.9259 , 



d\ = -9.07738 , d% = -1.42836 , d\ = 10.7003 , 



c\ = 60.9607 , 4 = 40.9998 , 4 = -0.348161 , 
c| = 3.21593, c% = 5.28888, 



4 = -0.715392 
4 = -475.465 , 
c b e = -347. 271 , 



4 = - 



490.982 , 
5.13592, 



f 2 = -1325.19, 



4 = -9.00634, 

4 = 47.1982, 
/£ = -219.224 , 



4 = 12.2908, c b 5 = -113.125, 
4 = 886.503 , 4 = -25.4831, 
4 = 91.1767, 4 = -297.002, 
/f = -283.955 , A fc = 736.209 , 
fl = 634.499 , fl = -25.8203 , 



fg = 301.478 , fl = -904.161 , / 6 Q = -271.966 , fl = 827.319 



4 = -0.0309341 , 

4 1 = -12.4700, 
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